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The chemical composition of ultra-high energy cosmic rays is a key question in particle astro- 
physics. The measured composition, inferred from the elongation rates of cosmic ray showers, looks 
in general very different from the initial source composition: resonant photo-disintegration in the 
cosmic radiation background proceeds rapidly at the highest energies and the initial composition 
quickly becomes lighter during propagation. For a statistical analysis of continuously improving 
cosmic ray data it is desirable to know the secondary spectra as precisely as possible. Here, we 
discuss exact analytic solutions of the evolution equation of ultra-high energy cosmic ray nuclei. 
We introduce a diagrammatic formalism that leads to a systematic analytic expansion of the exact 
solution in terms of second order effects of the propagation. We show how the first order corrections 
of this expansion can improve the predictions of secondary spectra in a semi-analytical treatment. 

PACS numbers: 96.50.S-, 98.70.Sa, 13.85.Tp 



I. INTRODUCTION 

The mass composition of ultra-high energy (UHE) cosmic rays (CRs) remains an open question in astrophysics. 
The average mass number (A) per energy can be inferred directly in CR observatories by the measurement of the 
elongation rate distribution of CR showers [T]. Presently, the observational situation is ambiguous despite strong 
experimental efforts over the years. Recent findings of the Pierre Auger collaboration [2] indicate a transition of UHE 
CRs within the energy range 10 18 eV to 4 x 10 19 eV from a light (presumably proton-dominated) spectrum towards 
a heavier composition [3]. In contrast, the HiRes collaboration 4 finds a mass composition compatible with that of 
a proton-dominated spectrum [5]. 

Various features in the CR spectrum can also provide indirect evidence for the origin and composition of UHE CRs. 
The ankle - a hardening of the spectrum at 3 x 10 18 eV - could be formed naturally by the superposition of two 
power-law fluxes and serves as a candidate of the transition between galactic and extra-galactic cosmic rays [51 [7] . It 
has also been advocated that this feature could be well reproduced by a proton-dominated power-law spectrum, where 
the ankle is formed as a dip in the spectrum from the energy loss of protons via Bethe-Heitler pair production (HE]. In 
this case extra-galactic protons could already start to dominate the spectrum beyond the 2nd knee which corresponds 
to a slight softening of the spectrum at 5 x 10 17 eV. 

Proton-dominance beyond the ankle is ultimately limited by the Greisen-Zatspin-Kuzmin (GZK) cutoff [101 111] 
due to resonant photo-pion production in the cosmic microwave background (CMB). In fact, a suppression of the CR 
spectrum at the expected energy of about 5 x 10 19 eV has been detected by the Pierre Auger and HiRes collaborations 
at a statistically significant level [H [T2] and is consistent with a proton dominance at these energies. However, this 
feature could also originate from photo-disintegration of UHE CR nuclei in the cosmic background radiation, or from 
an in situ energy cut-off of the injection spectrum of UHE CR. To summarize, the interpretation of these experimental 
findings is as yet inconclusive and even controversial. 

Simple theoretical arguments, however, can motivate a significant contribution of primary nuclei at energies beyond 
10 18 eV. For the efficient acceleration of primary CRs to these extreme energies a particle should be confined mag- 
netically in a suitable astrophysical environment. Since the particle's Larmor radius is proportional to its rigidity, 
i.e. its energy per charge, we expect that the maximal energy -E max of UHE CRs to scale with the charge number 
Z of a (fully ionized) nucleus. The acceleration of heavy nuclei like iron (Z = 26) can hence proceed up to larger 
energies and alleviates the fundamental limitations of cosmic accelerators to account for the observed spectrum of 
UHE CRs [13]. 

Analytic 1 descriptions of UHE CR propagation provide an easily accessible means of exploring both proton and 
heavy nuclei source scenarios. With the most recent results of the Auger collaboration indicating that the composition 



1 We use the term "analytic" here to denote explicit analytic solutions in closed form following [14] . There also exist many implicit 
analytic solutions for the spectra of UHE CRs that require an algorithmic treatment, e.g. [8l 1151 [16] ■ 
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continues to become heavier at energies above 3 x 10 19 eV, the heavy UHE CR flux component may well arrive from 
very local cosmological regions. To facilitate the future exploration of nearby UHE CR nuclei source distributions, 
we here develop further the analytic description of UHE CR nuclei propagation put forward in Ref. [13]. These 
developments take into account subdominant energy losses, ensuring that this description provides an accurate means 
of obtaining the UHE CR flux over the full energy range covered by the cosmic ray observatories. 

We will start in section [H] with a short recapitulation of the compact evolution equations of UHE CR nuclei in 
the limit of a spatially homogeneous distribution of isotropic CR sources. We derive an exact analytic solution in 
section III that includes continuous energy losses and multi-nucleon loss transitions between nuclei. In section IV we 
introduce a perturbative expansions of the exact analytic solution that provides a convenient practical framework for 
next-to-leading order corrections of the solution given in Ref. [14] . We finally conclude in section [V] 



II. PROPAGATION OF COSMIC RAY NUCLEI 



For a spatially homogeneous distribution of cosmic sources, emitting UHE particles of type i, the co-moving number 
density Yi is governed by a set of (Boltzmann) continuity equations of the form: 

Yt = d E {HEYi) + dEibiYi) - T\ ot Y t + j d£y 7,i^ + A , (1) 

together with the Friedman-Lemaitre equations describing the cosmic expansion rate H(z) as a function of red-shift 
z? The first and second terms on the r.h.s. describe, respectively, red-shift and other continuous energy losses (CEL) 
with rate b = AE/dt. The third and fourth terms describe more general interactions involving particle losses (i — > 
anything) with total interaction rate r' ot , and particle generation of the form j — > i with differential interaction rate 
7ij. The last term on the r.h.s., Ci, corresponds to the emission rate of CRs of type i per co-moving volume. 

The two main reactions of UHE CR nuclei during their cosmic evolution are photo-disintegration [T7H20] and Bethe- 
Heitler pair production [ST] with the cosmic background radiation (CBR). The former process is dominated by the 
giant dipole resonance (GDR) with main branches A —> (A — 1) + N and A —> (A — 2) + 2N where N indicates a proton 
or neutron [T71 - HU] . The GDR peak in the rest frame of the nucleus lies at at about 20 MeV for one-nucleon emission, 
corresponding to Eq DR ~ A x 2 x e~J v x 10 10 GeV in the cosmic frame with photon energies e = e me v meV. At 
energies below 10 MeV there exist typically a number of discrete excitation levels that can become significant for low 
mass nuclei. Above 30 MeV, where the photon wavelength becomes comparable or smaller than the size of the nucleus, 
the photon interacts via substructures of the nucleus. Out of these the interaction with quasi-deuterons is typically 
most dominant and forms a plateau of the cross section up to the photo-pion production threshold at ~ 145 MeV. 
Bethe-Heitler pair production can be treated as a continuous energy loss process with rate b^iz, E) — Z 2 b p (z, E /A), 
where b p is the energy loss rate of protons [2JJ. The (differential) photo-disintegration rate Ta->-b(E) ('Ja^b(E, E')) 
is discussed in more detail in Appendix [A} 

The evolution of the spectra proceeds very rapidly on cosmic time scales and the flux of secondary nuclei, J, looks 
generally quite different from the initial injection spectrum, Jj n j . The reaction network of nuclei depend in general on a 
large number of stable or long-lived isotopes. If the life-time of an isotope is much shorter than its photo-disintegration 
rate it can be effectively replaced by its long-lived decay products in the network ([I]). Typically, neutron-rich isotopes 
/3-decay to a stable or long-lived nucleus with the same mass number. In most cases there is only one stable nucleus 
per mass number below 56 Fe with the exception of the pairs 54 Cr/ 54 Fe, 46 Ca/ 46 Ti, 40 Ar/ 40 Ca and 36 S/ 36 Ar (see 
Fig. [9]) . We follow here the approach of Puget, Stecker and Bredekamp (PSB) [TH] and consider only a single nucleus 
per mass number A in the decay chain of primary iron 56 Fe. This PSB-chain of nuclei linked by one-nucleon losses is 
indicated as a red arrow in Fig. [9] 

As described earlier, CR nuclei that undergo rapid photo-disintegration with CMB photons carry a Lorentz factor 
of about 7 = 2 x 10 10 . We can only strictly neglect long-lived secondary isotopes from the reaction network if the 
nucleus lifetime in the cosmic frame, 7T, is much smaller than the inverse photo-disintegration rate, which is of the 



2 This is given by H 2 (z) = H 2 [Cl m (l + z) 3 + normalized to its value today of Ho ~ 70 kms — 1 Mpc -1 , in the usual "concordance 
model" dominated by a cosmological constant with Qa ~ 0.7 and a (cold) matter component, f2 m ~ 0.3 [13| . The time-dependence of 
the red-shift can be expressed via dz = — dt (1 + z)H. 
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order of (4/A) Mpc. This corresponds to nucleon life-times of less than a few minutes. Figure [9] shows also isotopes 
below 56 Fe with life-time larger than about one minute in addition to the nuclei of the PSB-chain. In general, there is 
a large number of isotopes that are sufficiently long-lived in the cosmic frame to take part in the photo-disintegration 
process. Fortunately, a large degeneracy of intermediate isotopes with equal mass number affects only very heavy 
nuclei. The photo-disintegration of these degenerate nuclei, dominated by collective excitations of nucleons like the 
GDR, mostly depend on the mass number A. The fluxes calculated for nuclei in the PSB-chain are expected to give 
a good representation of the total flux per mass number. Note that most of the analytic formulae that we are going 
to introduce in the following can be easily generalized to the case of the full reaction network including all isotopes. 

Note, that the Boltzmann equations ([I]) do not take into account the deflection of charged CR nuclei during 
their propagation through inter-galactic and galactic magnetic fields. The strength of inter-galactic magnetic fields 
is limited to the range 10~ 16 G - 10~ 9 G [531 H3] and suggested to be of 0(1CT 12 )G by simulations of large-scale 
structure formation |24j . In fact, if synchrotron radiation during propagation is negligible and the source distribution 
is homogenous, Eq. ([I]) provides a good approximation of the spectral evolution even for CRs having small rigidity 
which suffer large deflections |25j . However, magnetic inhomogeneities on small scales will suppress the spectrum of 
CRs with Larmor radius < where id is the characteristic distance between sources. It has been shown that for 
particularly strong inter-galactic magnetic fields of strength ~ 1 nG and coherence length of ~ 1 Mpc, the diffusive 
propagation of CR protons will start to affect the spectrum below about 10 9 GeV if l& ~ 50 Mpc [26] . Depending on 
the diffusion regime, this can suppress the proton flux at 10 8 GeV by a factor of 3 to 100. Due to the dependence 
£l oc 1/Z we expect that for heavy nuclei diffusive propagation can in principle remain important up to the ankle. 
The results of this paper are based on solutions of Eqs. ([I]) and assume that the contribution of inter-galactic or 
galactic magnetic fields can be neglected for the calculation of the UHE CR spectrum. 



III. ANALYTIC SOLUTION 



The secondary nuclei produced via photo-disintegration carry approximately the same Lorentz factor as the initial 
nucleus and the differential interaction rate in Eqs. ([!]) can be approximated as ja-+b(E, E') ~ T a^b{E)5{E' — 
(B/A)E). It is hence convenient to express the energy of a nucleus with mass number A and red-shift z as A(l + z)E 
where E denotes the energy per nucleon. Introducing the CR density per co-moving volume and nucleon energy, 
Naj = AEi(l + z)AYa(z, (1 + z)AEi), and corresponding emission rates, Qa.i = A(l + z)AEi£(z,A(l + z)Ei) we 
can re- write Eqs. dip in the compact form 3 



where we define the rates: 



^ T(A,i)^f(B,i)NA,i + E ^{B,i)^(A,i)NB,i + Qa.i ■ 



B<A 



pCEL 
1 A.i 



(A,i)->(A,i-l) 



b A (z,A(l + z)Ej) 
A(l + z)AE l 
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(A,i)^(B,i) 



A^B 



(z,A(l + z)E i ) 



(2) 



(3) 



Hooper et al. [T3] discussed an analytical solution of Eqs. (|2[) for one-nucleon losses in the limit = and 

Qaa = 0. In fact, the solution of a more general interaction network with generalized interaction rates ^(A,i)-^(B,j) 
of the form ^ can be written 



N A ,i(t) 



E E(nw]E 

j>i,B>A c 



1=1 



k=l 



d.t'Q B ,j{t')e 



-(t-t')rl ot 



n 



ptot _ ptot 

p=l(^fc) C P 



(4) 



where we sum over all possible production chains c = (ci, . . . , c„ c ) with intermediate nuclei of mass number C in the 
energy bin k - denoted by the doublet c$ = (C,k) - and fixed endpoints c\ — {B,j) and c„ c = (A, i). The partial 
width T C[ ^ C[+1 includes nucleon-disintegration (r^ ^^^ ^) as well as CEL (T^A,i)^(A,i-i))- A proof of Eq. Q is 
given in Appendix |B| 



3 This form of the differential equation holds for nuclei heavier than beryllium. We can easily compensate for the process 9 Be — > 4 He + 
4 He + n of the PSB chain (see AppendixjAl by re-defining N' A i = 7V4 j/2 for A = 2, 3, 4 and N' A . = Na,x f° r other nuclei. 
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FIG. 1: A possible transition chain c between an initial configuration (blue dot) and a final configuration (magenta dot) 
including one-nucleon losses (red arrows), two-nucleon losses (red dotted arrows) and continuous energy loss (green arrows). 
For the exact analytic solution Q all possible transition chains of this type are taken into account. 



We can visualize the production chains c diagrammatically as paths along the configuration grid of nuclei, as shown 
in Fig.[l] A horizontal link corresponds to a CEL transition whereas a vertical link denotes photo-disintegration. The 
color coding in Fig. [T] indicates the type of the transition q — > q+i - green for CEL, red for one-nucleon losses and 
red-dotted for two-nucleon losses 4 . It is convenient to use this graphical representation as a short-hand notation for 
the terms of Eq. Q. To see this, we can write Eq. Q in the form 

oo 

NaM= ( dt ' Yl G{A,i,B,r,t-t')[Q Bd {t') + 5{t')N Bd {Q)]> (5) 

{ j>i,B>A 

where we define a Green's function G(A, i, B,j; At) = &(At) J2 C G{c; At) as a sum over the contribution per path, 

(Tic — 1 \ Tic n c -t 

n r ™ E" ff ' r n (6) 
i=i / fc=i P =i( 7 tfc) c p ° k 

Each term G(c; At) in the previous equation corresponds to a production chain on the configuration grid. We will 
use this graphical representation later for a perturbative expansions of Eq. Q . 

The interaction rates T are not constant as the Universe expands. For example, the photo-disintegration rate 
with the CMB photons scales with red-shift as Ya{z,E) = (1 + z) 3 Ta(0, (1 + z)E), which follows from the adiabatic 
expansion of the CMB. Also, the nucleus emission rates La are not in general constant with time. A standard 
approach approximates the scaling with red-shift as a simple power-law over a finite red-shift distance, e.g. 

C A (z, E) = Q(z - z min )e(z max -z)(l + z) n C A (0, E) . (7) 

We can account for the red-shift dependence of V and Q by summing Eqs. Q over sufficiently small red-shift intervals, 
in which these quantities can be regarded as constant. Typically, intervals of Az ~ 0.01 are sufficient for this approach. 

Though the expression ^ is an exact analytical solution of the system of differential equations ^ , its calculation 
involves a large number of possible production chains and becomes numerically inefficient for large configuration grids. 5 
For instance, for one-nucleon and two-nucleon losses the number of possible chains Faa between nuclei with mass 
number A and B = A + A A can be derived iteratively from the identity Faa+2 — Faa + Faa+i with Fq = F\ = 1, 
which we recognize as the sequence of Fibonacci numbers. Hence, the total sum over different chains and primary 
nuclei in expression Q involves Fq + Fi + . . . + Fff—l = -Fjv+i ~ 1 number of terms, which is a number that scales 
exponentially with N. Hence, considering all transitions via one-nucleon and two-nucleon losses between, say, proton 
(A = 1) and iron (^4 = 56) becomes numerically very expensive even without considering transitions via CEL. 



4 We will use later on "generalized" chains, where the transition c; — > cj +1 is not necessarily equal to r ci _n, i+1 . 

5 In general, the numerical evaluation of expression (jSJ requires a high computational precision. We use the publicly available multiple 
precision libraries GMP 27 and MPFR 28 for this purpose. 
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FIG. 2: A graphic representation of the NLO paths contributing in the first correction for AA = 3 (see Eq. (|14||). The 
black arrows indicate transitions between configurations with total transition rate F* 
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and the dotted red arrows 



two-nucleon transition rates T A u respectively. Note that these types of graphs contribute with opposite sign in Eq. (h4h 



We show in the following that the exact expression Q can be well approximated by the dominant production chain 
through one-nucleon losses. Corrections via two-nucleon losses and CEL can be treated perturbatively. As means of 
a comparative check, we obtain results using our analytic description, assuming a source injection spectrum of the 
form 

J inj cx £^e- £ / £ — . (8) 
These analytic results are compared against those obtained numerically through a Runge-Kutta method [25] . 



IV. PERTURB ATIVE APPROACH 



The dominant contribution to the nucleon transitions in the CRB comes form one-nucleon losses with transition 
rate T^. In the following we focus on perturbative corrections to this dominant decay route from the contributions 
of two-nucleon losses and CEL with transition rates and T^ L , respectively. 



A. Two-Nucleon Losses 



We start with perturbative corrections from two-nucleon losses and assume, for the moment, that T^* = 
and = 0. For a perturbative expansion it is convenient to rewrite Eq. (B2) as 



p2N 
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(9) 



with 



ABC 



En 



i=i 



ptot 
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ptot 
1 C,i 
ptot 
1 DA 



(10) 



We can define a perturbative expansion of Eq. ( 10 1 in terms of sub-dominant branching ratios of two-nucleon pro- 
duction, r^/r^?'. The leading order (LO) contribution, J^^bc = 1> reproduces the approximation of Ref. [13] . The 
next-to-leading order (NLO) contribution can be written as 
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p2N 
1 D.i 
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(11) 
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FIG. 3: The solution Q at leading-order (LO) and up to next-to-leading order (NLO) (Eq. (14|) for one-nucleon (IN) and 
two-nucleon (2N) losses. To aid the comparison between the results, we ignore the evolution of the nucleon emission rates and 
interaction rates with red-shift and sum over red-shift steps Az = 0.01. We compare the LO and NLO analytic results to a 
numerical solution via a Runge-Kutta method |29| . 



We can most easily visualize these terms by a perturbative expansion of the nucleon densities, 
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where the terms N A [ are solutions to the set of differential equations, 
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(13) 



For the moment, we assume that the total width T^* is the sum of one-nucleon and two-nucleon losses. This, however, 
can be generalized to the total photo-disintegration rate for general nucleon losses (see section IV B). As an initial 



condition we define iV^ n ](0) = for n > and JV^(0) = N A i (0). Note, that with this initial condition the expansion 

(12 1 becomes finite and hence converges trivially. Each term corresponds, by construction, to the rt-th order 
correction of T . We can write the NLO correction explicitly as 

t 
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(14) 

Inserting the LO solution in Eq. (14 1 and following similar algebraic steps as in Appendix |b| one can identify .M 1 ' as 



the difference of contributions form paths {{A, i), . . . , (£?, i)) with length B — A + 1 and B — A, respectively, with the 
single insertion of a two-nucleon loss step into the decay chain. This is displayed diagrammatically in Fig. [2] for the 
case AA = 3. 
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Note, that we can also express Eqs. ([14]) as a matrix equation of the form, 

A^j(Ai) ~ J2 (XABAAt)N B 40)+y AB 4At)Q B ^) . (15) 



B>A 



The matrices X{At) and y(At) are in general only slowly changing with the red-shift scaling of the background 
radiation. It is hence possible to improve the NLO results by introducing sufficiently small time intervals Az and 
apply Eq. ((TBI repeatedly. 



We show the LO and NLO results of our approach in comparison to a numerical solution via a Runge-Kutta method 
in Fig. [3] For simplicity, we assume that CEL is absent and that source terms and interaction rates are constant 
throughout the integration domain < z < 1. The NLO contributions are shown for two cases. In the case "Az = 1" 



we calculate the NLO contribution directly by Eq. (14). The case "Az = 0.01" shows the improvement of the NLO 



contribution by a repeated application of Eq. (151 for the corresponding time interval - 100 times in this case. In 



most cases, the LO approximation is already satisfactory [k 



B. General Photo-Disintegration Losses 



For high mass nuclei (A > 40) of the PSB-chain one-nucleon and two-nucleon losses constitute more than 90% of the 
total photo-disintegration rate as can be seen in the Table [TJ However, for low mass nuclei the emission of a particles 
(as well as deuterons (D) and tritons (T)) can become important. As in the previous section, we can organize these 



sub- leading contributions via the expansion (12). For instance, the additional contribution from a particle loss can 
be introduced at NLO (n > 0) as 

J\r( n ) _ ptot Ar(n) , -ptot \r(n) , r 2N ivrt™- 1 ) p2N aj(«-1) , ro *r(™-2) r a J\r(«- 2 ) fl a\ 

iV A,i - ~ L A.i^A.i + 1 A+l,i JV A+l,i + 1 A+2,i^A+2A ~ 1 A+l,i iV A+l,t + 1 A+4,i ly A+i,i ~ 1 A+l,i JV A+l,i ' l 10 -' 

where we now have to include a emission in the definition of the total rate, T^' = + + T A ^ The treatment 
of these additional photo-disintegration channels is completely analogous to the case of two-nucleon losses. 



C. Continuous Energy Losses 



We next consider the contribution of CEL to the solution Q. In this case we have to include all possible paths in 
Eq. ( 4j that allow for both, variation of energy and mass number as the one shown in the left panel of Fig. [I] Similar to 
the discussion of two-nucleon losses, the number of possible paths becomes very large. For the remainder of this section 
we consider only one-nucleon photo-disintegration losses together with CEL and, hence, T^* = + r^ L - Despite 
this simplification there are still (AA + Ai)l / (AA)l / (Ai)l different paths in total between the two configurations 
(A, i) and (A + A A, i + Ai). This becomes computationally very expensive for long production chains, as we already 
observed for the introduction of two-nucleon losses. 

We can account for CEL transitions as effective source terms in the differential equations This turns out to be 



an efficient way for determining the resulting spectra. As before, we can use the perturbative expansion ( 12 ) of the 

(n) . ... 

nucleon densities, where the terms N A [ are now solutions to the set of differential equations, 

= ~^t<\ + r^iX+M + Q*,i . ( 17 ) 
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iV A,i - _i A,i ISI A,i + 1 A+l,i iV A+M + A,i+l^AA+l ~ 1 A+l,i iV A+l,i {U > U) . 

Here, the total width 1^°* = + is now for the sum of one-nucleon and CEL, though in general it would 

receive contributions from all exclusive channels. Note that with the initial condition -/V^ ?(0) = for n > the 



expansion (12 1 of N Ai is finite if a finite set of energy bins and nuclei is considered, A < A 



specifically, the expansion of N A ,% only includes non-zero terms for n < (A max + z max ) — (A + i). 
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FIG. 4: A graphic representation of the NLO paths contributing in the first correction N A . 
black arrows indicate transitions between configurations with total transition rate Ta.\ = 
transitions with CEL rate r^ L , respectively. Note that these types of graphs contribute with opposite sign in Eq. ( 18 \ 
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The first term Nj£\ in the expansion of Naa is our familiar solution for the one-nucleon loss case (|9|) where the 



"A,i 

partial width is replaced by the total width. The second term TV^ can be evaluated explicitly by an insertion of N^°> , 

t 



EE n rt ^ 



B>AC=A \D=A+1 



n 



i 



ptot ptot 

D=A(^C) D - i C - i 



(18) 

After some algebraic manipulations one can identify A^ 1 ) as the difference of contributions form paths 
((A, «),..., (B, i + 1)) and ((A, i), . . . , (B, i)) with the single insertion of CEL step into the decay chain. This is 
indicated diagrammatically in Fig. [4] for the case A A = 3. 

Note that the NLO correction for CEL only introduces transitions between the energy bins i and i + 1. Hence, 
the NLO solution ( 18 1 can not be considered as a small correction to the full solution if the contribution from CEL 
becomes large, Atr^^ > 1. However, in analogy to the case of two-nucleon losses we can write the NLO contribution 
as a matrix equation 



N ( 2}{At) ~ i^AB ,i (At) Ns,i (0) + y A B,i{M)Q B< i + V A B,i( At ) N B,i+l(0) + W AB ,i(M)QB,i+l) 



(19) 



B>A 



If we consider sufficiently small time intervals At such that Atr^ L « 1 we can approximate the exact solution 



by a repeated application of Eq. (19). The transition matrices X(At), y(At), V(At) and W(Af) are only slowly 
changing with the scaling of the background radiation. It is hence only necessary to re-evaluate these matrices on 
large time-scales; typically Az ~ 0.01 is sufficient for the propagation of heavy nuclei. Thus, results obtained by the 
application of this procedure should be considered semi-analytic. 

Figure [5] shows the results of the LO and NLO energy flux spectra compared with results obtained using a Runge- 
Kutta method. For simplicity, we again consider constant source terms and interaction rates and assume that two- 
nucleon losses are absent. The repeated application of Eq. (JT9J) reproduces the numerical solution well. For heavy 
nuclei (and hence "short" transitions from primary iron) or large energies E/A > 10 10 GeV the LO contribution is 
already an excellent approximation. 



D. Secondary Proton and Helium Spectra 



Finally, we discuss an expansion of the spectrum of primary and secondary protons 6 and helium. This case is 
slightly different from the propagation of heavy nuclei, since there are additional contributions from the channels 



We do not distinguish between protons and neutrons in the following, assuming a prompt decay of secondary neutrons. 
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FIG. 5: Comparison of the terms in expression Q up to next-to-leading order (NLO) with the numerical solution via a Runge- 
Kutta method including one-nucleon (IN) and continuous energy losses (CEL). To aid the comparison between the results, we 
ignore the evolution of the nucleon emission rates and interaction rates with red-shift and sum over red-shift steps Az — 0.01. 




FIG. 6: The full NLO correction for two-nucleon and continuous energy losses in comparison with the numerical solution. For 
this result we took into account source evolution and red-shift effects, assuming that the nucleon emission rates scale as (1 + z) 3 
and sum over red-shift steps Az — 0.01. 
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FIG. 7: An example of an n-th order production path contributing to N^f including the effective nucleon production rate F°bj 
(blue arrow). The box in the top right corners indicate the complete sum over all possible n-th order contributions "J^ 



Mr 



to 



the production chain of the nucleon (B,j) with B > 2 and j > i. The red arrows indicate the CEL contribution for protons, 
including energy loss by Bethe-Heitler pair production and photo-pion production. 



(7,N), (7,2N), (7, a), (7,Na) and (7,2a). Secondary nucleon production follows the differential equation 
where the effective nucleon production rate from transitions (A,i) (l,i) is defined as 

nos,N _ t^in 1 2r 2N + r Nct + 5 " plN 1 * " p2N ' 1 * ^NO-I 



r 



AA 



AA 



AA 



r 



AA 



), 



(20) 



(21) 



with = 1 if ^4 = £? and zero otherwise 7 . Note, that the last term in (21) is assumed absent in the PSB-chain 



Photo-hadronic interactions of the protons can be determined using the Monte Carlo Package SOPHIA 30 . Here, we 
approximate photo-pion interactions of the protons as a continuous energy loss process in addition to Bethc-Hcitlcr 
pair production. The differential equation (20) is of the same form as Eq. (2) and we can hence write its exact solution 
in the form Q. 

With the expansion ^ we can write the set of evolution equations as 



r CEL m O , r CEL aHO , peff at n 

1 1A iV l,i + 1 M+l iV M+l + 1 AA 1SJ AA + ' 



(22) 



A>2 



N 



(n) 



A>2 



tCEL »t(«) 



eft ^(n) 



(n > 0) 



In contrast to the case of nuclei, we cannot treat CEL of the protons as a second order effect. Nevertheless, with 
the set of differential equations (p2k and the boundary condition N^(0) — iVi^O) and n[ " (0) = for n > the 
expansion (pi) is finite since the expansion of Nai is finite. Explicitly, we can write the n-th order contribution as 



EE IK! L )E 

B>2 j>i \k=i+l 



k—i 



n 



1 



pCEL _ pCEL 



(23) 



Again, these contributions to the proton spectra can be expressed via diagrams indicated in Fig. [7j By definition, the 
term N^f depend on all possible n-th order production chains of intermediate nuclei (B,j), that are indicated as the 
boxes in the top right corner of the diagrams. 

Similarly, the emission of a particles in the channels (7,0?), (7,2a) and (7,Na) rate can be described by the differential 
equation 



N 4 



r^iiv 4 . 



r?f L iV4,i 



(Tft + I%)Nt,i + 52 rff^i + Qaa ■ 



(24) 



A>2 



For N' A . = N A ,i/2 for A = 2, 3, 4 we re-define Ff / = 2Ff v 
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FIG. 8: Left: The average mass number from a pure-iron _E _2 -flux with £ max = 10 22 eV and source evolution parameters 
n = 3 and z max = 1. We show the full numerical solution in comparison with the LO and NLO analytic equation including 
the exclusive channels (7, N), (7, 2A), (7,0), (7, Not) and (7,2a) as well as CEL. Right: The total energy flux for the same 
parameters plotted against recent Auger measurements [2]. The LO results are in good agreement with those shown in Fig. 4 
of [M]. 



with an effective production rate 

T Tf = T A.i + T A." + 2F A° + ^12 r !2 :l + $A9Fg°; (+S A 8^8,i) ' ( 25 ) 



Again, the last term in (251 is absent in the PSB-chain considered in our calculation. In principle, we can treat these 
contributions analogously to the case of the protons. However, the relative contribution from a particle emission is 
only small if we consider heavy primary nuclei like 56 Fe and can be neglected in this case. 

The sum over diagrams of the type shown in Fig. [7] involve a large number of intermediate configurations (£> , j) 
and the calculation can become time-consuming. For a more efficient calculation of the proton spectra we can utilize 
the total conservation of nucleons per energy bin within sufficiently small time-steps with AtTff L <C 1. In this case 
the flux can be well approximated as 

N hl (At) ~ JV M (0) + AtQi,i + At [r^A^+^O) - T^N^Q)] 1 N aM + ^Qa,% ~ N A;i (At)} . (26) 

A>2 

With this approximation, and using the NLO contribution of the exclusive channels (7, AT), (7, 2 A 7 "), (7, a), (7, Na) 
and (7, 2a) as well as CEL for the spectra of nuclei, we show in the left panel of Fig. [8] the average mass number (A) 
in comparison with the analytic result. The right panel of Fig. [8] shows the total energy flux of nuclei for the NLO 
analytic solution compared to the numerical result. For these results, time steps of Az = 1CP 4 have been used in order 
for the proton contribution to the total flux to be calculated with the necessary accuracy. The LO approximation is 
already in excellent approximation to the data of CR observatories considering the large systematic and statistical 
uncertainties of the CR spectra and the average mass composition. All spectral features of the quantities and their 
overall scale are well reproduced by the LO contribution. Improvements to the LO result, however, are made by the 
NLO contributions, whose results leave only a very mild discrepancy with the Runge-Kutta results at energies below 
!()" ' GeV. 



V. CONCLUSIONS 



In this work we have developed further an analytic solution for the fluxes of UHE CR nuclei from extragalactic 
sources. We have shown that in most cases the spectra are well approximated by the analytic solution already given in 
Ref. [14], which dealt with the dominant energy loss channel of single nucleon transitions between nuclei. We have here 
expanded on this approach through the introduction of NLO corrections from two-nucleon and CEL. The introduction 
of these terms was shown to further improve the accuracy of the analytic description. In order for these results to 
take into account the slow variation of interaction and emission rates with red-shift as well as CEL we incorporated 
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our result into a semi-analytic framework. The semi-analytic results obtained were found to be in excellent agreement 
with results obtained through a purely numerical Runge-Kutta approach. 

The prospects of determining the nature of extragalactic UHE CRs and their sources in the near future are promising. 
Ongoing direct hybrid measurements of UHE CRs by the Auger collaboration continue, with the opportunity now 
existing for an independent verification of these results by other hybrid experiments such as the Telescope Array |31j . 
These measurements allow the possibility for a coherent picture of the UHE CR flux, composition, and arrival direction 
anisotropy to emerge. Present and ongoing indirect measurements of the secondary particles produced by UHE CRs 
during their acceleration and propagation are also capable of constraining the UHE CR composition and their sources. 
For instance, the simultaneous emission of neutrinos arising from proton-proton and/or proton-photon interactions in 
extra-galactic protons sources can serve as a test of low energy crossover scenarios [8j of extra-galactic protons [16j 122] . 
Photo-pion interactions by extra-galactic protons in the CMB, i.e. the process responsible for their GZK-cutoff, give 
rise to a flux of cosmogenic neutrinos [3"3"H3"5] and photons [3BJ . The accompanying output into secondary electrons 
and positrons, in particular from Bethe-Heitler pair production, feeds into electromagnetic cascades in the cosmic 
background radiation and intergalactic magnetic fields [37]. This leads to the accumulation of 7-rays at GeV-TeV 
energies. The observed extra-galactic diffuse 7-ray flux thus provides a constraint on the total energy injected into 
such cascades over the Universe's entire history |38| . 

The methods provided in this paper offer a general tool with which theoretical results may be easily obtained and 
compared to both these direct and indirect UHE CR measurements. As example cases, the application of the general 
methods developed here to proton propagation provide the opportunity to further develop the method applied in 
[39) . Secondly, an analytic determination of the photon fraction produced through UHE CR nuclei propagation is 
anticipated to also be obtainable using this treatment. Through the simplicity of our approach and the speed with 
which it may be implemented, our analytic method is anticipated to be of great benefit as a tool for future UHE CR 
investigations. 
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Appendix A: Photo-Disintegration of Nuclei 

The most general evolution of primary and secondary nuclei in the CRB includes all possible photo-disintegration 
transitions between nuclides (A, Z) competing with the decay of unstable nuclides. For simplicity, we follow the work 
of Puget, Stecker and Bredekamp (PSB) [r7TfT9l and consider only one stable isotope per mass number A in the decay 
chain of 56 Fe as already explained in section |ll| This "PSB-chain" is listed in Tableland sketched in Fig. [9} 

Table |TJ shows the relative contribution of inclusive channels to the total photo-disintegration rate calculated for the 
nuclei of the PSB-chain. We use the reaction code TALYS [20] to evaluate the cross sections for nuclei with 10 < A < 56 
and assume an E~ 2 power-law flux of CR nuclei. At CR energies E < 10 12 GeV and large mass numbers A > 20 
photo-disintegration in the CRB can be well approximated by one-nucleon and two-nucleon losses between elements 
of the PSB-chain via exclusive processes (j,p), (7, n), (7, 2p), (7, 2n) and (7,2m). For the cross sections of light nuclei 
with mass numbers A = 2, 3, 4 and 9 we use the parameterization of Ref. [40] . 

At lower mass numbers, A < 20, additional channels involving a particle emission can become as significant as the 
sume of one-nucleon and two-nucleon losses. Table|TJalso shows the relative importance of the exclusive channels (7, a), 
(7,71a), (7, pa) and (7,2a) to the total photo-disintegration budget. Resonant photo-nuclear interactions play only a 
minor role in the propagation of the nuclei for the energies of interest. We follow the approach outlined in Ref. [40] 
and approximate the total interaction by the isospin averaged rate as T J 4 7 (z, E) ~ ATnj(z, E / A). We also assume 
that the participating nucleon is removed from the nucleus and regard this as a contribution to one-nucleon losses. 

The angle-averaged interaction rate appearing in Eq. (JTJ) is then defined as 

1 

T A ^ B {z,E) = ^ J dcosfl J de(l- Pcose)n 1 {z,e)a A ^ B {e), (Al) 
-1 
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FIG. 9: The Puget-Stecker-Bredekamp-chain 18 along with stable and long-lived (r > 1 min) nuclei below 56 Fe. 



where n 7 (z,e) is the energy distribution of isotropic background photons at red-shift z and e' = 67(1 — j3 cos 9) the 
photon's energy in the rest frame of the nucleus. For our calculation we use the cosmic microwave background and the 
infra-red/optical background form Ref. @T]. To a good approximation the decay products of the photo-disintegration 
interaction inherit the large boost-factor of the initial nucleus and hence in the process A —> B + (A — B) the nucleus 
with mass number B has an energy E' — (B/A)E. We can hence approximate the differential cross section as 

~/a^b(E, E') ~ T A ^ B {E)8((B/A)E - E') (A2) 

in the following. This has the correct normalization since Ta-+b(E) = j dE'jA^siE, E'). 

In general, the interaction rates r a^b{z, E) scale with red-shift according to the red-shift evolution of the radiation 
background. In the case of the CMB with adiabatically scaling, n 7 (z, e) = (1 + z) 2 n 7 (0, e/(l + z)), we can derive the 
simple relation 

T A ^ B {z,Ei) = (l + zfT A ^ B (0 1 (l + z)E). (A3) 

For the case of the infra-red/optical background [41] we assume a red-shift scaling following the star formation rate 
as described in Ref. [T5]. However, since the cascades of UHE CR nuclei develop locally, the red-shift dependence of 
the interaction rates is only of minor importance. 
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TABLE I: The nuclei of the Puget-Stecker-Bredekamp-chain [18] and the relative contribution of inclusive channels to the total 
photo-disintegration cross section in the CMB calculated by TALYS [2U]. We assume and E~ 2 spectrum of the nuclei and 
integrate over nucleon energies 10 17 eV < E/A < 10 21 eV. Channels with contribution less than 1% are omitted in the table. 

Appendix B: Proof of Equation ^ 

We will proof Eq. Q by induction. First note, that we can rewrite Eq.Q as 

n 

= -r\ ot Ni + Fj^iNj + Qi , (Bl) 

3=1+1 

with i — 1, ... ,n with Tj^ t = for i < j and r' ot 7^ rJ- ot for i ^ j. In the following we will refer to the indices i 
as knots and the pairs (i, j) with ^ as links. A chain of length n c is defined as an ascending sequence of n c 
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knots, ci < C2 < . . . < c na , that are mutually connected by links. 



We want to show that the most general solution of Eq. (|Bl| is of the form 

/n c — 1 \ n c 



j<i c \ i=l 



fe=l 



^■(0)e- tr ^ + 



d*'Q i (t / ) e - (t - t ' )r " 



TT 1 

J. 1 ptot _ ptot 

o=l(/fc) c p Cfc 



where the sum is over all possible chains c with c\ = j and c„ c = i. 
Induction start: n — 1. This case has the solution 

t 

N x (t) = JVi(0)e- tr f t + / dt'Q^e-^'^ . 



(B2) 



This is of the form (B2 ), since the only chain is the trivial one of length n c — 1 with c\ = 1. 



Induction step: n n + 1. The differential equations of AT, with 1 < i < n are of the form (Bl I and we can hence 
use the solution (B2). The differential equation for N n+1 is 

n 

N n+1 = -r^iJV„+i + £ r m ^ n+1 iV m + Q„+i , (B3) 



We can write the general solution of this differential equation as: 



iV n+ i(t) 



JV n+1 (0)e- tr "+i + / dt'Q n+1 {t')e-^'^ 



P n 

+ J dt' e -(*- t ') r "+i £ r m ^ n+1 JV m (t') 



(B4) 



The first term of the previous equation corresponds to the first term (i = j = n + 1) in the sum of Eq.(B2 1. Inserting 
the solutions (B2) in the integrand yields after integration by parts: 



N n+1 (t) 



N n+1 (0)e- tT ^ + / dt>Q n+1 (t>)e-^i°U 



(B5) 



( n r^^+i ) r m^n+i^ i p tot _p tot I r iHi _ r i ( „ 



m— 1 j — 1 c \ / — 1 



^(0)e~ tr 'r + / dt'Q^t^e-^'^ 



Nj (0)e- tT ^ + / dt'Q^e-^ 1 '^ 



The chains c in the previous sums have end-points c\ = j and c ric — m. Now, every chain c in the system with n 
knots and endpoint c„ c = m corresponds unambiguously to a chain c' in the system with n + 1 knots with = q for 
i < n c i — 1 and c' n r = n + 1. Hence, the double-sum in Eq.(B5 1 over end-points m < n and chains c can be expressed 
as a single sum over chains c' with d-y = j and d n l = n + 1. We arrive at the form: 



N n+1 (t) 



iV„ +1 (0)e- tr "+ 1+ / dt'Q n+1 (t')e-^'^ 



(B6) 



n /n c (-l \ 

EE n r <-< +1 ) £ 

j = l c' \ 2 = 1 



fc=l 



iVj(0)e °i + / dt'Q 3 (t')e 



fr 1 

X X ptot ptot 

p=l(/fc) C P c 'fc 



/II,/ -1 



j = l c' \ 1 = 1 



Ni(0)e-^ + / dt'Qj^e-V-*'^ 



n a i — l n c i 

e n 

k=i P =i(^k) 



i 



ptot ptot 
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As a final step we use the identity: 8 



n c / — 1 n c / | n c i — X 

5z n ptot _ ptot _ _ rr 



r 



tot ptot 



(B7) 



n+l 



to combine the last two terms in Eq. ( B6 ) and arrive at the form ( B2 ) 
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